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Abstract. We present a massive equilibrium simulation of the three-dimensional 
Ising spin glass at low temperatures. The Janus special-purpose computer has allowed 
us to equilibrate, using parallel tempering, L = 32 lattices down to T w 0.64Tc. 
We demonstrate the relevance of equilibrium finite-size simulations to understand 
experimental non-equilibrium spin glasses in the thermodynamical limit by establishing 
a time-length dictionary. We conclude that non-equilibrium experiments performed 
on a time scale of one hour can be matched with equilibrium results on L w 110 
lattices. A detailed investigation of the probability distribution functions of the spin 
and link overlap, as well as of their correlation functions, shows that Replica Symmetry 
Breaking is the appropriate theoretical framework for the physically relevant length 
scales. Besides, we improve over existing methodologies to ensure equilibration in 
parallel tempering simulations. 
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1. Introduction 

Spin Glasses (SG) are disordered magnetic alloys that are generally regarded as 
particularly convenient model systems for the study of glassy behaviour [H [2] . Indeed, 
ideas originating in the SG context have been fruitful in the study of structural glasses, 
optimisation in computer science, quantum information, econophysics, etc. 

A distinctive feature of SG is that, below their glass temperature, they remain out 
of equilibrium even if they are left to relax under constant experimental conditions for 
days or weeks. In spite of this, the equilibrium properties of their low-temperature phase 
is believed to control their non-equilibrium behaviour. Indeed, both theory [3l H] and 
experiment [5] agree in that the sluggish dynamics is due to a thermodynamic phase 
transition at a critical temperature T^., that separates the paramagnetic phase from a 
low-temperature one where the spins freeze according to extremely complex, essentially 
unpredictable, ordering patterns. Furthermore, it has been now established that an 
accurate knowledge of the thermodynamic equilibrium properties would allow us to 
predict in detail many relevant features of their non-equilibrium relaxation [HI [7]. 

There is an already 30-year-old theoretical controversy regarding the defining 
properties of the SG phase. On the one hand, the Replica Symmetry Breaking (RSB) 
theory that stems from Parisi's solution of the SG in the mean field approximation [8l[9]. 
A system well described by the RSB is in a critical state for all T < Tc, where the 
surfaces of the magnetic domains are space filling. On the other hand, the droplet 
theory [HI [121 113] views the SG phase as a disguised ferromagnet. It provides the 
solution of SG models as computed in the Migdal-Kadanoff approximation [H]. We 
refer the reader to section [2] for the detailed predictions of the RSB and droplet theories 
for the different physical observables in the SG phase. The predictions of the somewhat 
intermediate TNT theory [15], [16] are discussed also in section [2l 

Numerical simulations are the main tool that theoretical physicists have to make 
progress in the understanding of the SG phase in D = 3 systems. Basically without 
exceptions, numerical work in Z) = 3 is best described by RSB theory (see [9] for a review, 
refs. [T71 [IHl [19] for recent work and refs. [151 CSl 120] for some somewhat dissenting 
views). Yet, numerical investigations have received as well severe criticism. It has 
been claimed that basically all simulations doable to date are contaminated by critical 
effects [21]. One would need to simulate still larger systems at still lower temperatures, 
in order to observe the asymptotic behaviour corresponding to large enough systems. 

Here we present the results of a large-scale simulation performed on Janus [221 123] , 
a special-purpose computer designed for the simulation of SG. For this particular task, 
Janus outperforms standard computers by several orders of magnitude. We have devoted 
(the equivalent of) 200 days of the full Janus computer to an equilibrium, parallel- 
tempering simulation of the Ising SG in D = 3. We have been able to thermalise lattices 
of size L = 32 down to temperatures T ^ 0.64Tc. This is not only a world record, but 
provides as well an unprecedented glimpse on the low temperature SG phase. 

Our main objectives here have been (see section [2] for definitions): 
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• To perform a precision comparison of equilibrium and non-equilibrium spatial 
correlation functions. It turns out that a time-length dictionary exists, which 
relates with amazing accuracy our previous results at finite times [211 125] (on non- 
equilibrium infinite systems) with equilibrium finite lattice sizes. The unavoidable 
conclusion is that experimental SG are in the dynamical non-equilibrium regimes 
that correspond to equilibrium results on lattices L ~ 110. There is no doubt 
that at these length scales, the appropriate effective theory is RSB, irrespectively 
of which of the competing theories is correct for much larger L. 

• To perform a study of the probability density function (pdf) of the spin overlap, 
and to extrapolate important quantities to the thermodynamic limit. So doing, 
we will gather important information about the correlation length in the spin-glass 
phase. 

• To provide a detailed study of the link overlap. 

• Last, but not least, to obtain a large set of configurations, fireproof thermalised, 
which will serve as a starting point for more sophisticated studies (such as 
investigation of ultrametricity, or temperature chaos). In particular, a detailed 
study of the spatial correlation functions will appear elsewhere [26] . 

The layout of the rest of this paper is as follows. In section [2] we briefly recall 
the definition of the Edwards- Anderson model. In particular, in section [2l2] we describe 
the observables considered and discuss the scaling behaviour predicted for them by the 
different theoretical scenarios. In section [HI we describe our simulations and address 
the crucial problem of ensuring thermal equilibrium. We have found it most useful 
to study the random walk in temperature space performed in our parallel-tempering 
simulations (section \3.3\\ . In particular, our thermalisation checks significantly expand 
the methodology introduced in [27]. At this point, we are ready to study in section H] 
the pdf of the spin overlap. In particular, in section I4l3] we determine through finite size 
effects a correlation length in the spin-glass phase. We focus on the spatial correlation 
functions in section [5l finding (section [6]) crystal-clear indications of the relevance of our 
equilibrium investigations to the non- equilibrium experimental work. The properties of 
the link overlap are addressed in section [71 Our conclusions are presented in section [HI 
Technical details are provided in two appendices. 

2. Model, Observables, Theoretical expectations 

We divide this section in five paragraphs. In section [2m we describe our model. The spin 
overlap and related quantities are defined in section 12.21 We discuss spatial correlation 
functions in section 12.31 Their non-equilibrium counterparts are recalled in section 12.41 
We address the link overlap in section 12.51 Even though most of this section consists of 
results and definitions well known in the spin-glass community, we consider it convenient 
as a quick reference. We also introduce some specific (and sometimes new or seldom 
used) physical quantities for this paper. 
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2.1. The model 

We consider the D = 3 Edwards- Anderson model [28l[29]. Our dynamical variables are 
Ising spins Sx = ±l, which are placed on the nodes, x, of a. cubic lattice of linear size 
L, containing V = sites, and with periodic boundary conditions. Their interaction 
is restricted to lattice nearest neighbours and is given by the Hamiltonian: 



Note that the couplings Jx,y in the Hamiltonian are themselves stochastic variables: 
they take the values ±1 with 50% probability. The coupling constants attached to 
different lattice links are statistically independent. The physical motivation for working 
with a random Hamiltonian is modelling the effects of impurities in a magnetic alloy. 

We shall consider the quenched approximation: in the time scale relevant to the 
spin dynamics, the impurities can be regarded as static. Hence, we will not allow for any 
back- reaction of the spins over the coupling constants. A given realisation of the {Jx,y} 
(a sample, from now on), will be fixed from the start and considered non-dynamical [1]. 

A random Hamiltonian implies a double averaging procedure. For any observable O 
(an arbitrary function of the spins and the coupling constants) , we shall first compute the 
thermal average (O) using the Boltzmann weight at temperature T for the Hamiltonian 
(II]). The average over the coupling constants distribution, (O) , is only taken afterwards. 
We will refer sometimes to the second averaging, (■ ■ ■), as disorder average. 

The reader will notice that the disorder average induces a non-dynamical gauge 
symmetry [30]. Let us choose a random sign per site = ±1 . Hence, the energy ([T]) is 
invariant under the transformation 



Since the gauge-transformed couplings ex^yJx,y are just as probable as the original 
ones, the quenched mean value of {0{{sx])) is identical to that of its gauge average 
X]{e^=±i} {0{{exSx})) /'2.^° , which typically is an uninteresting constant value. We show 
in section 12.21 how to overcome this problem. 

We remark as well that the Hamiltonian ([1]) also has a global Z2 symmetry (if all 
spins are simultaneously reversed Sx —Sx the energy is unchanged), corresponding to 
time-reversal symmetry. This symmetry gets spontaneously broken in three dimensions 
upon lowering the temperature at the SG transition at = 1.109(10) [311 132] . 

2.2. The spin overlap 

We need observables that remain invariant under the transformation (|2]). The 
Hamiltonian ([T]) provides, of course, a first example. To make further progress we 
consider real replicas {sx^}, {sx^}, copies of the system that evolve under the same set 




(1) 



Jr. 




(2) 
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of couplings {Jx,y} but are otherwise statistically uncorrelatedI|| 
Using them we form the overlap field: 

= s^^^s^:^ , (3) 

which is obviously invariant under (|2]). 

The Edwards- Anderson order parameter, the spin overlap, is the spatial average of 
the overlap field: 



X 



In particular, it yields the (non-connected) spin-glass susceptibility 



XNc(T) = V{q^) , (5) 

that diverges at with the critical exponent 7. For all T < Tc, one expects 
Xnc = ^iy) ■ We shall also consider the Binder ratio 



B{T) = ^, (6) 

In particular, for all T > T^, the fluctuations of q are expected to be Gaussian in the 
large-L limit, hence lim^^-^^oo = 3, (T > T^). Its behaviour in the low-temperature 
phase is controversial. For a disguised ferromagnet picture one expects B to approach 
1 in the limit of large lattices. On the other hand, for an RSB system one expects 
1 < i? < 3 in the SG phase (T < Tc). We recall also that one may consider as well the 
overlap computed in small boxes, in order to avoid the effect of the interphases (physical 
results are equivalent to those obtained with the standard overlap [33]). 

A great deal of attention will be devoted to the probability density function (pdf) 
of the overlap 



X 

Note that, in a finite system, the pdf is not smooth, but composed of + 1 Dirac deltas 
at g = — 1, • • • , "^1^5 1- Here, we have solved this problem by a convolution of the 

comb-like pdf d?]) with a Gaussian of width 1/W, Qvix) = J ^exp[— "K^] : 



P{q = c) = P dq' P{q')Gy{c-q') = (^Gy{c-^Y.'i-))- 



In this way, we basically add the contribution of 0{yV) microscopic values of g, 
belonging to an interval of width ~ l/\/V [3l]. Note, however, that eqs. ( 15)161) are 
computed out of moments of -P(g), rather than of P{q)- 

The Edwards- Anderson order parameter ^ea vanishes for all T > Tc. Below Tc, in 
a droplet system, P{q) collapses in the large-L limit in a pair of Dirac delta functions 

I For the thermal average of any observable depending on a single spin configuration, 0({si^^}), we 
have (0({si'^}))' = (0({si'^})0({si')})). 
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of equal weight, centred at g = ±gEA- In an RSB system, P{q) contains as well a pair 
of delta functions at q-^A, but it also has a continuous piece, non- vanishing for every q 
such that — gEA < q < qEA- This is the origin of the differences in the predictions that 
both theories make for B in the low-temperature phase. 

We will find it useful to consider as well conditional expectation values at fixed q. 
Let O be an arbitrary function of the spins. We define its conditional expectation 



Of course, one may easily recover standard expectation values from E(0|g): 



(0)= / dgP(g)E(0|g). (10) 

J — oo 

Strictly speaking, the integration limits should be ±oo. However, truncating the integral 
to — 1 < g < 1, the error is exponentially small in L^^'^ (yet, for L = 8 and 12 we had to 
extend the limits beyond ±1). 

We can also define the conditional variances as 

Var(0|g = c) = E{0'^\q = c) - E(0|g = c)^ (11) 

where we have the identity 



(02) - (O) = / dg P(g) Var(0|g) + (E(0|g) - {0)Y . (12) 

J — oo 

2.3. Spatial correlation functions 
The overlap correlation function is 

X 

Ci^{r) decays to zero for large r only for T > T^. Thus we have considered as well 
conditional correlation functions, recall eq. (|9]): 



C4(r|g) = E 




g . (14) 



Eq. ( ITOl) allows us to recover C^^r) from C4(r|g). 

The two main theoretical pictures for the SG phase, the droplet and RSB pictures, 
dramatically differ on their predictions for C4(r|g). Let us discuss them in detail: 

• In the RSB picture, the connected correlation functions tend to zero at large r. For 
all g G [— gEA7Q'EA] we expect the asymptotic behaviour 
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where the dots stand for scahng corrections, subleading in the hmit of large r. The 
exponent 6{q) in eq. (fT5|) has been computed for D larger than the upper critical 
dimension = 6: [351 EE] 

e{q = 0) =D-A, (16) 

^(0<|g|<gEA) = /^-3, (17) 
d{\q\ = qEA) =D-2. (18) 

These mean-field results for 6{q) become inconsistent for D < 4 [the correlations 
should decrease for large r, implying 6{q) > 0, recall eq. (fTSl) ]. An expansion in 
e = 6 — D suggests that 6{q) will renormalise [37]. Note as well that, at least for 
large D, 6{q) is discontinuous at g = 0. However, we remark that there are no 
compelling theoretical arguments supporting the discontinuity of 9{q) in D = 3. 
Indeed, recent numerical studies found no evidence for it [191 126] . We finally recall 
a non-equilibrium computation [25] yielding in D = 3l§| 

^(g = 0) = 0.38(2). (19) 

• Quite the opposite to the RSB case, in a system well described by a droplet model 
and for |g| < ^ea, ^4(^1?) does not tend to for large r (we are referring, of course, 
to the regime 1 <^ r <^ L). In fact, spin configurations with \q\ < ^ea are spatially 
heterogeneous mixtures of the two pure phases. One should find bubbles or slabs of 
linear size ~ L of one of the two phases, say q = +qEA, surrounded by a matrix of 
the complementary state (see e.g. [38l [39] ). It follows that 

C4{r\q) = qlJr/r{r/L) , if |g| < gEA and K r < L , (20) 

{fr/r{x) is a direction- dependent scaling function with /r/r(0) = 1). Indeed, 
the probability that two spins at fixed distance r belong to domains of opposite 
orientation is proportional to r/L in the large-L limit. On the other hand, precisely 
at |g| = qEA but only there, droplet theory predicts that the connected correlation 
function vanishes for asymptotically large r. The same behaviour of eq. (|T5l) 
was predicted [11]. The exponent 0{qEA) is identical to the scaling exponent of 
the coupling strength, denoted as 6* or y in the literature, and has a value of 
^(gEA) ~ 0.2 [n]. 

2.4- Non- equilibrium correlation functions 

Let us recall that non-equilibrium counterparts exist of g and C4(r|g). We shall not be 
computing them here, but we will compare previous computations with our equilibrium 
results. Hence, we briefly recall the definitions [25]. One considers pairs of times tw and 

§ Wc may mention as well three conjectures: 6{0) = {D — 2 + rj)/2 [37] (that from the results 
in [32], yields 9{0) = 0.313(5)), 6{0) ~ \/v {v is the exponent that rules finite size effects at (/ea) 
and 0(0) + l/j) = 9 {qEA)- There is also an exact scaling relation 0{q-EA) = 2/i> [26] , 
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t + 1„, with t, t„ > 0, after a sudden quench from a fully disordered state to the working 
temperature T. The analogous of the spin overlap is 

C{t, ^w) = :^ XI Mt^)Sx{t + t^)) . (21) 

X 

The non-equilibrium spatial correlation function is 

C2+2(r; t, tw) = {SxiU)Sx{t + t^)Sx+r (t + tw)) (22) 

X 

At fixed tw, C{t,tvi) monotonically decreases from C = 1 at t = 0, to C = at t — )■ oo. 
Hence, one may consider C, rather than t, as an independent variable. We will compare 
the non-equilibrium C2+2{r]t,t^), computed in very large lattices [211 [25], with our 
equilibrium results for C4 (r|g = C{t,t^)). To do so, we shall need to relate the finite 
time tw (on very large lattices) with the finite size L. As we shall see in section [6|, the 
correspondence between the non-equilibrium and the equilibrium correlation functions 
is amazingly accurate. 



2.5. The link overlap 

The link overlap is defined a^ 

Qiink = ^ Yl ^^^v ■ 

||a;-y|| = l 

It is a more sensitive quantity than the spin overlap to the differences between a system 
described by droplet theory or an RSB system [40]. Since it is invariant under time- 
reversal symmetry (the global reversal of every spin in either of our two real replicas 
s^x — > —s^x) its expectation value is non- vanishing, even in a finite system at high 
temperatures. Its pdf can be defined as we did with the spin overlap, recall eqs. ( ITIIHI) . 
In fact, it has been proposed that the link overlap (rather than the spin overlap) should 
be considered as the fundamental quantity to describe the spin-glass phase below the 
upper critical dimension ^\\ [T7] . There are both physical and mathematical reasons for 
this: 

• On the physical side, Qimk provides an estimate of the volume of the domains' 
surfaces. Indeed, consider two configurations of the overlap field ([3]) differing only 
in that a domain of size ~ L has flipped. This will result in a large change of the 
spin overlap, q. Yet, the only changing contribution to Qunk is that of the lattice 
links crossed by the domain's surface. In a droplet theory, where the surface-to- 
volume ratio of the domains vanishes in the large-L limit, one does not expect any q 
variation of the conditional expectation E(Qiink|?), not even in the \q\ < qEA region. 
Hence, the pdf for Qunk is expected to collapse to a single-valued delta function in 
the large-L limit. The intermediate TNT picture coincides with the droplet theory 
in this respect. For an RSB system, the domains' surfaces are space filling. Hence, 



II Clearly, (Qii„k) = ^(1, 0, 0). 
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when q suffers a variation of order 1, the variation of Qunk will be of order 1, too. 
Accordingly, a non-trivial pdf is expected for Qunk? in the limit of large systems. 

On the mathematical side, theorems have been proven for the link overlap |l2l |33l 
m] , valid for three-dimensional systems, which are the exact correlate of mean-field 
results for the spin overlap^ Specifically, the replica equivalence property holds 
for the link overlap in three dimensional systems. Replica equivalence [45l |46] is a 
property of the Parisi matrix which yields an infinite hierarchy of identities relating 
linear combinations of moments of Qimk in the large-L limit. A specific example 
that we shall be using here is 



lim (Qiink)^ = lim 

L—^oo L— >oo 



1 



- (Qlink) + - (Qunk) 



(24) 



(at finite L, the equality is not expected to hold). This is just a particular case of 
the family of identities valid for all /c, s = 0, 1, 2, ... 

1 



lim (QLk)(^?fi„k)= lim 



2 (Qfink) ('^llnk) + 2 (Qvmk) 



(25) 



(replica equivalence implies infinitely many relations such as this). It is amusing 
that the mathematical proof for the three-dimensional theorem does not use Parisi 
matrices, relying instead on stochastic stability. Let us stress that ultrametricity 
implies replica equivalence, but the converse statement (i.e. replica equivalence 
implies ultrametricity) does not hold, in generalQ 

The distinction between spin overlap and link overlap seems somewhat artificial 
from the point of view of the mean-field approximation. In fact, in the Sherrington- 
Kirkpatrick model one easily shows that Qunk = q"^- For finite-connectivity mean- field 
models, non-equilibrium numerical computations yield Qimk = OQ'^ + b [IE] (a and b are 
numerical constants). In D = 3 there are also clear indications that fixing the spin- 
overlap fixes as well the link overlap: the conditional variance Var((5iink|9), eq. ( fTTl) . 
tends to zero for large lattices, see [17] and figure [161 below. Furthermore, in a TNT 
or droplet system, the derivative dE{Qii^]^\q)/dq'^ should vanish in the large-L limit for 
all |g| < gEA (since there is a single valid value for Qunk, there can be no dependency 
left). Numerical simulations, both in equilibrium [171 [H] and out of equilibrium [211 HQ], 
find so far a non-vanishing derivative that nevertheless decreases for larger L. The 
extrapolation to L = oo is still an open issue, see section 17.11 

We wish to emphasise that Qunk unveils that the spin-glass phase is a critical state 
where minimal perturbations can produce enormous changes. In fact, let us couple two 

% The mathematical proof known so far is vahd only for Gaussian-distributed couplings in eq.([T]). 
However, physical intuition strongly suggests that the theorems are valid in more general cases such as 
our bimodal couplings. 

+ For the sake of completeness, let us recall that replica and overlap equivalence, combined, imply 
ultrametricity . In addition, replica equivalence and the Ansatz of a generic ultrametricity implies 
ultrametricity just as in the SK model |47| . Finally, we point out that replica equivalence is tantamount 
to stochastic stability and a self-avcrageness property. 
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otherwise independent copies of the system through Qunk, 



n 



- TeVQ 



link 



In a system described by droplet theory, one expects the hnk susceptibihty 



Xlink 



5(<5link) 



de 



V 



(Qfink) - {Q 



link 



(26) 



(27) 



to remain finite in the large-L hmit, for all T < Tc (precisely at Tc, a critical divergence 
might arise). Hence (Qiink)^ = {Q\ink)^=o + exiink + • • • in a droplet or TNT system. 

On the other hand, in the mean-field approximation, one finds for RSB systems a 
discontinuity with e [50] : 

= E((5iink|g = ^ea) + a+^/e+..., (28) 
= E(gunkk = 0)- a_v^+... . (29) 



(Qlink) 



e>0 



(Qlink)e<o 

Actually, the mean-field computation was carried out for the spin overlap, yet, in 
mean-field models, Qwnk is essentially g^, hence we can borrow their result. We should 
emphasise that the situation is even more critical than for standard first-order phase 
transitions: Xiink(e) diverges when e — (just as if the specific heat of liquid water 
approaching its boiling temperature showed a divergence!). 

Below the upper critical dimension, there has been very little investigation of xiink 
(see, however, ref. [IQ]). In fact, eq. (l24|) has interesting implications in this respect. Let 
us rewrite it in the equivalent form 



lim 

L— >oo 



(QLk) - (^?link)^ 



- lim 

3 L— >oo 



{Q 



2 ) 
link/ 



(Qlink) 



(30) 



In an RSB system, the right-hand side of eq. ( I30l) is positive (since Qunk may take values 
on a finite interval). Yet, eq. (127|) . the Ihs of ( 130|) is nothing but the large-L limit of 
X\ink/L^- Hence, RSB implies xunk ~ L^, as expected for first-order phase transitions 
(see e.g. [5Tj). 

We note that for droplet, or TNT systems, eq. fl30|) is merely an empty = | x 
statement, just as for RSB systems in their paramagnetic phase. Hence we have found 
of interest to study the dimensionless ratio 



R 



link 



{Q 



2 ) 
link/ 



(Qlink)' 



('^Hnk) ~ (Qlink) 



(31) 



eq. fl30|) implies that, for an RSB system on its large-L limit, -Runk = § for all T < T^,. 
For a droplet or TNT system any value < -Runk < 1 is acceptable. In fact, the 
high-temperature expansion for the D = 3 EA model tells us that, in the large-L limit, 
i?u„k = l-0(T-2). 

We finally recall that the Chayes et al. bound [521 153] may seem to imply that xunk 
can diverge at most as L"^/^, rather than required by RSB. The way out of the 

paradox is a little technical O 

* Imagine generalising model ([T]) in the following sense: the coupling is J^y = +1 with probability p 
(and Jxy = — 1 with probability 1 — p). so that our model is just the particular instance p = 0.5. One 
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3. Numerical methods 

We describe liere our numerical simulations. We describe the simulation organisation 
on Janus in section 13.11 We explain our choice of parameters for the parallel tempering 
simulation in section 13.21 An absolutely crucial issue is that of thermalisation criteria, 
section 13.31 We largely extend here the methods first introduced in ref. [27], which 
allows us to distribute on a rational basis the computational resources according to the 
difficulty in thermalising each particular sample. At variance with ref. |27j, which was 
restricted to the critical region, we are here probing the deep spin-glass phase, hence 
more demanding criteria need to be met. The statistical data analysis is described in 
section 13.41 Finally, in section 13.51 we describe some more traditional thermalisation 
tests. 

3.1. The Janus computer 

Our Monte Carlo simulations have been carried out on the Janus special-purpose 
machine. Information about Janus' hardware as well as some details of low-level 
programming can be found in [22], |23l [Ml- Janus is built out of 256 computing cores 
(Virtex-4 LX200 FPGAs) arranged on 16 boards. With the code used for this paper, 
each core updates 3 x 10^*^ spins per second with a heat bath algorithm. The 16 FPGAs 
on a board communicate with a host PC via a 17th on-board control FPGA. 

The controlling PC generates the couplings {Jxy}, initialises the Janus random 
number generators, and provides as well the starting spin configurations. All the 
required data is transmitted to the FPGAs (one FPGA per real replica) that carry 
out both the Heat Bath (HB) updating and the Parallel Tempering (PT) temperature 
exchange. Due to the special architecture of Janus, the PT step is not costless, as we 
previously need to compute the total energy for each temperature. We thus equilibrate 
the computational cost of both updates by performing several HB sweeps before a PT 
temperature swap is attempted. Fortunately, selecting a modest number of HB sweeps 
per PT update hardly affects the efficiency. After a suitable number of PT cycles, spin 
configurations of all replicas are copied to PC memory to take measurements. The 
measurement process on the PC is easily parallelised with the next simulation block in 
Janus so that the PC is always ready for the next reading. 

During the simulation, we store on disk information about the PT dynamics 
(temperature random walk and acceptance rates), configuration energies, and 
measurements related to the overlap and link overlap fields. We also store full spin 
configurations every several measurement steps (usually a hundred) to be later used for 
offline measurements (see section 13.41) or as a checkpoint for continuing the simulation 

may follow ref. [52] to show that d{Q\ink) /dp diverges at most as L^/^. However, the eritieal value of 
e would still be e = for p in a finite range around p = 0.5 (this is the crucial point: in the standard 
argument [52[ 153] one would require that the critical value of e vary when p moves away from p = 0.5). 
Hence, the rate of divergence of p-derivatives does not convey information on the rate of divergence of 
e-derivatives. 
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if needed. 

In a few specific cases (namely one L = 24 sample and four L = 32 samples) the 
time required to fulfil our thermalisation criteria was exceedingly long, more than six 
months. For these samples we have accelerated the simulation by increasing the level 
of parallelism. We have used a special low-level code that transfers the PT procedure 
to the control FPGA. This has allowed us to distribute the set of temperatures along 
several FPGAs on a board, speeding up the simulation accordingly. 

For the smaller lattices {L < 12) we substitute the communication with Janus 
by a call to a simulation routine in the PC. Although these simulation are much less 
demanding, we go down to very small temperatures. As a consequence, the total cost 
is not negligible and we have used a PC cluster to complete the simulations. 

3.2. Choosing parameters for Parallel Tempering 

The key point in a parallel-tempering [55l |56] simulation consists in ensuring that each 
configuration spends enough time at high temperatures so that its memory can be 
erased. Since we intend to study the physics of the Edwards-Anderson spin glass at 
very low temperatures, our simulations are necessarily very long. Because of this, we do 
not need to reach temperatures as high as those used in critical point studies. We can 
perform a quantitative analysis using the known behaviour of the heat-bath dynamics 
above the critical point. 

Following [57], the equilibrium autocorrelation time in the thermodynamic limit is 
taken from a power law to a critical divergence 

rHB(T)~(T-T,)-^^ (32) 

For instance, for the maximum temperature used in our largest lattice (L = 32) Ogielski 
found tub{T) ~ 10^ [57]. This is several orders of magnitude shorter than our shortest 
simulations (see tabled]). 

The choice of the minimum temperature was taken so that the whole simulation 
campaign took about 200 days of the whole Janus machine and so that T^, — Tmin ~ 
L"^/''. With 4000 samples for L = 16,24 and 1000 for L = 32, this resulted in 
^min = 0.479,0.625 and 0.703, respectively. Smaller lattices, L = 8, 12, were simulated 
on conventional computers. In all cases, we simulated four independent real replicas per 
sample. 

As to the other parallel-tempering parameters, namely the number and distribution 
of intermediate temperatures and the frequency of the parallel tempering updates, the 
choice is more arbitrary. We dedicated several weeks of the machine to test several 
combinations trying, if not to optimise our decision, at least to avoid clearly bad choices. 

Specifically, we varied the number of temperatures keeping the acceptance of 
the parallel-tempering update between 7% and 36%. This corresponds to an increase of 
roughly a factor of two in Nt. Noticing that the computational effort is proportional to 
Nt, we found that the efficency hardly changed, even for such a wide acceptance range. 
Eventually, we chose a compromise value of about 20% in the acceptance, resulting in 
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Table 1. Parameters of our parallel-tempering simulations. In all cases we have 
simulated four independent real replicas per sample. The Nt temperatures are 
uniformly distributed between Ty^in and T^ax (except for the runs of the first row, 
which have all the temperatures of the second one plus T = 0.150 and T = 0.340). In 
this table Na-^cs is the number of Monte Carlo Steps between measurements (one MCS 
consists of 10 heat- bath updates and 1 parallel-tempering update). The simulation 
length was adapted to the thermalisation time of each sample (see section 13.31) . The 
table shows the minimum, maximum and medium simulation times (iVue) for each 
lattice, in heat-bath steps. Lattice sizes L = 8, 12 were simulated on conventional 
PCs, while sizes L — 16, 24, 32 were simulated on Janus. Whenever we have two runs 
with different Tmin for the same L the sets of simulated samples are the same for both. 
The total spin updates for all lattice sizes sum 1.1 x 10^". 



L 


T ■ 

mm 


T 

max 


Nt 


'mes 


N^t 




N^B 




System 


8 
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the parameters quoted on table [H This both avoided unconventionally low acceptances 
and saved disk space. 

In contrast to conventional computers, Janus needs about as much time to do a 
parallel-tempering update than a heat-bath one. Therefore, while it is customary to 
perform both updates with the same frequency, after testing frequencies in the range 
1-100 we have chosen to do a parallel-tempering update each 10 heath-bath ones. In 
fact, even if the time to do a parallel-tempering step were negligible, we have checked 
that doing a single heat-bath between parallel temperings would produce a practically 
immeasurable gain. We note, finally, that this issue was investigated as well in ref. [58] 
(in that work clear conclusions were not reached, as far as the D = 3 Edwards- Anderson 
model at low temperatures and large L is concerned). 

3.3. Thermalisation criteria 

In order to optimise the amount of information one can obtain given a computational 
budget, the length of the simulations must be carefully selected. It is well known that 
sample-to-sample fluctuation is the main source of statistical error. Thus, we want to 
simulate each sample for the shortest time that ensures thermalisation. 

The most common robust thermalisation check consists in the determination of 
the autocorrelation times for physical observables [59]. However, in order for this 
determination to be precise one needs a much longer simulation than needed to 
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thermalise the system (e.g., while ten exponential autocorrelation times can be enough 
to thermalise the system, we need an at least ten times longer simulation to determine 
this autocorrelation time). Notice that this is not an issue in ordered systems, where 
one employs very long simulations in order to reduce statistical errors. 

The typical practical recipe to assess thermalisation for disordered systems consists 
in studying the time evolution of the disorder-averaged physical observables. In 
particular, the so-called log2-binning procedure uses the evolution of the time averages 
along the intervals = (2~("'+^^A'hb, 2~"A^hb]- The system is considered to be 
thermalised if the first few intervals are compatible. 

This procedure is not optimal, because the thermalisation time is wildly dependent 
on the sample. Thus, a simulation time long enough to thermalise the slowest samples 
will be excessive for most of the rest. Perhaps even more frightening, the average over 
samples may well hide that a few samples, the very worst ones, are still quite far from 
equilibrium. 

Fortunately the use of parallel tempering presents us with the possibility to use 
the dynamics of the temperature random walk to learn about the thermalisation scale 
for each sample. In fact, in order to ensure thermalisation each of the participating 
configurations must cover the whole temperature range. Here, expanding on a method 
first used in [27], we have promoted this idea to a fully quantitative and physically 
meaningful level. 

Let us consider the ordered set of A^^ temperatures {Ti, . . . , Tjy^} and let us suppose 
that Ti^-i < Tc < Ti^. In figure [1] — left we show an instance of the random walk of 
the temperature index, i{t) G {1,2,..., Nt}, performed by one of the A''^ copies of the 
system considered in the parallel tempering. The random walk is clearly not Markovian, 
as the system remembers for a long time that it belongs to the high (low) temperature 
phase. This effect is also demonstrated in figure [1] — right, where we plot the time spent 
over Tc as a function of the simulation time (mind the long plateaux). 

To make these arguments quantitative, we shall use the standard tools of correlated 
time series [511 EH]- We need a mapping defined on the 1, . . . , Nt range of temperature 
indices so that 



It is also convenient that / be monotonic. Because we have chosen the same number of 
temperatures above and below T^, a simple linear / is suitable, but the method works 
with any function fulfilling the above conditions. 

For each of the participating configurations we consider the time evolution it of the 




(33) 
(34) 




(35) 
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Figure 1. We plot (left panel) the temperature index of a fixed configuration of an 
L = 32 sample as a function of the number of HB sweeps. We plot one point every 
5 million HB sweeps. The critical temperature corresponds to ic = 17. This specific 
sample has Tcxp = 1-75 x 10^° HB sweeps. In the right panel, we show the time that 
all configuration of a same replica spends in the paramagnetic phase. 



temperature index. We define the equilibrium autocorrelation function as 

C{t) = -— (36) 

^HB - to - t ^tTo 

where to is long enough to ensure that the temperature random walk has reached a 
steady regime. Due to condition (l35l) . we avoid subtracting the squared mean value of 
/ in this definition. From the normalised C{t) = C(t)/C{0), see e.g. figure [21 we can 
define the integrated correlation times: 

1 ^ 

nnt = - + J2^it), (37) 

t=0 

where is a self-consistent window that avoids the divergence in the variance of Tint- 
The great advantage of these functions over the physical observables is that we can 

average over the A^^ configurations in the parallel temperingl^ 

This procedure works surprisingly well, not only giving reliable estimates of the 

integrated time but even providing the, more physical but notoriously difficult to 

measure, exponential autocorrelation time. Indeed, the correlation function admits an 

tt Even if these are not completely statistically independent, the averaged autocorrelation has a much 
smaller variance. In addition, the need to simulate several replicas provides independent determinations 
of C(t), which permits a further error reduction and an estimate of the statistical errors. 
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Figure 2. Autocorrelation functions for samples with Toxp of different orders of 
magnitude. We have plotted the range [0,6rexp]- We include the automatic double 
exponential fit, see [Appendix A[ In the last panel the fit fails due to the strong 
downwards fluctuation and our programme has chosen a restricted interval for a fit to 
a single exponential. In order to avoid cluttering the graphs, we have only plotted a 
few times (the actual correlation functions have many more points). The horizontal 
axis is in units of 10^ heat-bath updates. 



expansion on exponentially decaying modes 

Cit) = J2a, e-*/---% 5^ = 1. (38) 

i i 

In this representation, the exponential time Texp is the largest of the r^^p j lffl Barring 
symmetry considerations, this exponential time should be the same for all random 
variables in the simulation, including the physical observables. 

The relative sizes of the Ai, and hence Tint, depend to a certain extent on the 
particular choice of /. Notice, however, that criteria (l33l - [34|) select a family of functions 

ft The number of modes equals the dimension of the dynamical matrix of the Monte Carlo Markov 
process, which in our case is (NtI) x 2^^^. 
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Figure 3. Histogram of exponential autocorrelation times for our simulations of the 
i = 32 lattice (1000 samples). 
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Figure 4. Logarithm of the histogram of exponential autocorrelation times for our 
simulations of the L = 24 lattice (4000 samples). Mind the behaviour of the long-times 
tail. 



that hopefully reduce the amplitude of the irrelevant fast modes. In any case, Tcxp has 
a physical meaning independently of these somewhat arbitrary considerations. 

In practice, the simulations are too long (up to A'^hb ~ 10^^) to consider all the 
f{it) individually and we have to introduce some data binning, averaging over a large 
number of consecutive measurements. As it turns out, this is not a very limiting issue 
for two reasons. On the one hand, as long as these bins are much shorter than r, there 
is no real information loss. On the other hand, one can reconstruct any polynomial / 



Nature of the spin-glass phase at experimental length scales 



18 



up to degree k — in particular our linear / — by saving the sums of the first k powers 
of the if. 

Even after this binning, we have worked with time series with a length of up to 
several million, so in order to compute the autocorrelation we have used a Fast Fourier 
Transform algorithm |60] . 



The details of the chosen thermalisation protocol can be found in Appendix A 



We summarise by saying that our main thermalisation criterion is ensuring that 
A^'hb > 12rcxp (2rcxp are discarded and the remaining lOrcxp are used to measure and 
study C{t)). 

In figure |2] we plot several autocorrelation functions showing how the data quality 
allows for an exponential fit. We have chosen randomly 4 samples with very different 
exponential autocorrelation times: 6.5 x 10^, 8.8 x 10^, 1.5 x 10^ and 1.8 x 10^°. 

To summarise the distribution of the exponential autocorrelation times we have 
computed a histogram. Due to the large dispersion of these quantities we have chosen 
log2 Texp as a variable. In figure |3] we show the results for the two runs performed in 
L = 32 (see table [1]). Notice the dramatic increase of the Texp when decreasing the 
minimum temperature of the simulation. The smooth shape of the curves defined by 
the histogram is a further test of our procedure for determining autocorrelation times. 

In figure H] we plot the logarithm of the histogram in the L = 24 case to 
show the exponential behaviour of the long-times tail. This result gives confidence 
that rare events, with very large (logarithms of) autocorrelation times, are at least 
exponentially suppressed. We have not made efforts to measure with precision the 
small autocorrelation times as they are immaterial regarding thermalisation, which is 
ensured by the minimum number of iterations performed for all samples. 



3.4- Monte Carlo evaluation of observables 

We present now some technical details about our evaluation of mean values, functions 
of mean values and error estimation. 

Some of the observables considered in this work were obtained by means of an online 
analysis: the internal energy, the link overlap, powers of the spin overlap {q,q'^,q^), 
and Fourier transforms of the correlation function C^lr) for selected momenta. These 
quantities were computed as Monte Carlo time averages along the simulation. Note that 
the length of the simulation is sample dependent, something that would be a nuissance 
in a multispin coding simulation, but not in Janus were each sample is simulated 
independently. The disorder averaging followed the Monte Carlo one. Statistical errors 
were computed using a jackknife method over the samples, see for instance |51] . 

However, when designing the simulation, one cannot anticipate all quantities that 
would be interesting, or these can be too expensive to be computed in runtime. 
In particular, we did not compute the conditional correlation functions C4{r\q). 
Fortunately, an offline analysis of the stored configurations has allowed us to estimate 
them. We had to overcome a difficulty, though, namely the scarcity of stored 
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configurations. In fact, for the samples that were simulated only for the minimum 
simulation time, we had only A'conf ~ 100 configurations stored on disk (ranging from 
-^conf = 10 for L = 12 to Nconi = 200 in the case L = 32). We regard the second 
half (in a Monte Carlo time sense) of these configurations as fireproof thermalised. Yet, 
when forming the overlap field, eq. ([3]), one needs only that the two spin configurations, 
{sx^} and {sx^}, be thermalised and independent. Clearly enough, as long as the two 
configurations belong to different real replicas and belong to the second half of the 
Monte Carlo history they will be suitable. There is no need that the two configurations 
were obtained at the same Monte Carlo time (as it is done for the online analyses). 
Furthermore, the four real replicas offer us 6 pair combinations. Hence, we had at least 
6x (A''conf/2)^ ~ 10000 (60000 for L = 32) measurements to estimate the overlaps and the 
correlation functions. We used the Fast Fourier Transform to speed up the computation 
of the spatial correlations. For those samples that had more configurations (because their 
total simulation time exceeded A^^in); considered nevertheless A'conf/2 configurations 
evenly spaced along the full second half of the simulation. When some quantity, for 
instance the P{q), could be computed in either way, online or offline, we have compared 
them. The two ways turn out to be not only compatible, but also equivalent from the 
point of view of the statistical errors. As an example of this let us compute the following 
quantity: 

^link=Mj-(^'- (39) 

For L = 32, T = 0.703, the value of af^^^^ computed from online measurements of Qunk 
and Qfink is 

^^ikonline = 50.88(90). (40) 

We could now recompute this value from offline measurements of Qunk and Q^nk- 
Instead, we are going to use eq. (fT2l) . which involves the intermediate step of computing 
conditional expectation values and variances at fixed q and then integrating with P{q). 
This will serve as a test both of the offline measurements' precision and of our Gaussian 
convolution method for the definition of clustering quantities. The result is 

^^ikconf = 50.81(90), (41) 

The precision of crfi^]^ online and CTj^inkconf is the same and the difference less than 10% of 
the error bar, even though we only analysed 100 configurations per sample for the second 
one. Of course, both determinations are very highly correlated, so the uncertainty in 
their difference is actually much smaller than their individual errors. Computing the 
difference for each jackknife block we see that 

V[4nKconi - <k,onlino] = -0.065(79), (42) 

which is indeed compatible with zero. 

A subtle point regards non-linear functions of thermal mean values that are later 
on averaged over the disorder. In this work, the only instance is Xunk, see eq. ( 1271) . 
Care is needed to estimate such non-linear functions because a naive evaluation would 
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be biased, and the bias might be sizeable compared to the statistical errors ^T\. This 
problem does not arise in non-linear functions such as eq. (|9]), which are computed on 
observables only after the double averaging process over the thermal noise and over the 
samples. The problem and several solutions are discussed in Appendix B| (see also [32]). 

A final issue is the comparison of data computed in different system sizes at the same 
temperatures. Unfortunately the grids of temperatures that we used for the different L 
differ. Hence we have interpolated our data by means of a cubic spline. 



3.5. Thermalisation tests 

We will consider in this subsection thermalisation tests directly based on physically 
interesting quantities. 

We start with the traditional logj-binning procedure. We choose the Binder 
parameter for the overlap, see eq. ([6]), which is specially sensitive to rare events. In 
figure [5] we show the results for B{Tj^i^) for L = 32, considering only the first 4 x 10^ 
Heat Bath steps of each of our 1000 samples, as if all the simulations were A^^in heat- 
bath steps long (blue line). We could not affirm that even the last two bins were stable 
within errors. Things change dramatically if we consider Monte Carlo histories of a 
length proportional to the exponential autocorrelation time. Note that, thanks to our 
choice of A^^in table [Tj the simulation time for most samples has not increased. If 
we first rescale data according to the total simulation length (itself proportional to the 
autocorrelation time) and average for equal rescaled time, the log2-binning procedure 
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Figure 5. Evolution of the Binder parameter for L = 32, T = 0.703 using log2 
binning (0 = second half, 1 = second quarter, ...). The blue curve (circles) is the 
result of stopping at step 1 of our thermalisation protocol (i.e., all samples simulated 
for a fixed time of 4x 10^ heat-bath updates). The red curve (squares) is the result 
of completing all the steps, which implies an increase of roughly 150% in simulation 
time. 
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Figure 6. Binder ratio as a function of the temperature for L = 32. The good overlap 
between two different simulations (one of them in the much easier critical region) is a 
further thermalisation check. We use the same set of 1000 samples. 

gives 4 steps of stability within errors. That is to say: we obtain the Binder parameter 
without thermahsation bias just discarding 1/16 of the history (and taking up to 1/8). 
Regarding the Binder parameter our requirement of 12rexp is excessive. 

In retrospect (see figure [5]), shorter simulations would have produced indistinguish- 
able physical results for most observables. We do not regret our choices, however, as we 
plan to use these thermalised configurations in the future [62] for very delicate analyses 
(such as temperature chaos), which are much more sensitive to thermalisation effects. 

A different test can be performed by comparing the difficult low-temperature 
simulations of our largest lattice with simulations in the critical region of the same 
samples. A faulty thermalisation (for instance, a configuration remains trapped at 
low temperatures) could be observable as inconsistencies in the values of quantities in 
common temperatures. In figure [6] we show the Binder parameter as a function of 
temperature for the two simulations with L = 32 (see tabled]). The agreement between 
both simulations is excellent. 

A very different test on the statistical quality of our data is the comparison of the 
values of ximk obtained using the different possible estimators for (Qiink)^- We have an 
unbiased estimator if we use (5iink4R) see eq. ( IB. 61) . the linearly bias-corrected estimator 

<5iink,iinear ^q. (Ell), and the quadratically bias-corrected estimator Q{2k,quadratic in 
eq. flB.5|) . The different determinations are equal only if the total simulation time (in 
each sample) is much longer than the integrated autocorrelation time for Qunk- As we 
see in figure [7]-left, only computing xunk from the biased estimator [Qiink]2/2 results in 
a measurable bias. Once bias correction is taken into account, differences are only a 
fraction of the statistical error for each estimator. Nevertheless, the different statistical 
estimators are dramatically correlated. Hence, their difference might be significant. 
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Figure 7. Bias correction in the computation of xiink, eq. ([27| . On the left panel 
we plot the two-replica estimators xiink,2R as a function of the unbiased four-replica 
estimator xiink,4Rj eq. ()B.10|) . for all our temperatures in the L = 32 lattice. The 
two-replica estimators xiink,2R are computed with no bias correction, eq. (|B.7p . with 
linear corrections, eq. (|B.8|) . and with quadratic corrections, eq. (|B.9p . The right 
panel displays, for the three two-replica estimators, their difference with the four- 
replica estimator in units of the statistical error for that difference, as a function of 
temperature. We show our data for L = 32 and L = 24. Note that the statistical error 
in the difference between two estimators is largely reduced (as compared to individual 
errors) due to dramatic data correlation. 



In figure [3-right we plot these differences for L = 24 and L = 32 as a function of 
temperature, in units of the statistical error for that difference. As we see, at the lowest 

(2) 

temperatures for L = 32, the bias for the estimate of xunk obtained from QiinkUncar 

(2) 

still measurable. Only the estimate from Qunk quadratic statistically compatible with 
the unbiased estimator. Since our data fully complies with our expectations, we consider 
the above analysis as a confirmation of our expectation S> Tint.Qunk) ''exp • 

We carefully avoided to make decisions during thermalisation based on the values 
of physical quantities. However, one could worry about the possibility of important 
statistical correlations between the temperature random walk and interesting quantities. 
Such correlation could originate some small biases that would be difficult to eliminate. 
Fortunately, we have not found any correlation of this type. In figure |8] we show the 
correlation between r^^p and two important quantities: probability of the overlap being 
small and the energy. 



4. The overlap probability density 

In this section we study the pdf of the spin overlap. This is a particularly interesting 
quantity because, as we saw in section |2l it has a qualitatively different behaviour in 
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Figure 8. Scatter plot of the exponential autocorrelation time [L ~ 32) versus the 
probability of the overlap being less than a small quantity (left) and the energy (right). 
We do not observe correlation between the thermalisation times and these physically 
relevant quantities. 

the droplet and RSB pictures of the spin-glass phase. 

We have plotted P{q) for T = 0.703 (the lowest for L = 32) and T = 0.625 (the 
lowest for L = 24) in figure |9l Notice that the convolution of the comb- like P{q), 
eq. (I7j), with the Gaussian function, eq. ([H]), has yielded a very smooth P{q). Initially, 
one would expect the peaks of this pdf to grow narrower and closer together as L 
increases, eventually becoming two Dirac deltas at ±gEA- The shift in position is clearly 
visible in the figures, but a more careful analysis is needed to confirm that the peaks 
are indeed getting sharper (section I4l3l) . In addition, the probability in the q = sector 
should either go to zero (droplet) or reach a stable non-zero value (RSB). Even if a 
visual inspection of figure [9] seems to favour the second scenario, we shall need a more 
quantitative analysis to draw conclusions. 




Figure 9. Overlap probability density function P{q), eq. ([5]), at T = 0.625 and 
T = 0.703. Notice that for the central sector oi q ^ the curves for the different 
system sizes quickly reach a plateau with P{q) > 0. 
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Figure 10. Overlap density distribution function at zero overlap as a function of 
temperature. We observe an enveloping curve with a linear behaviour, as expected in 
an RSB setting. 



In the remainder of this section we undertake such a quantitative characterisation 
of P{q) and, in particular, its thermodynamical hmit. To this end, we will study the 
evolution of P{q = 0) with T and L (section l4.ip : the extrapolation to infinite volume 
of the Binder cumulant (section I4.2p and finally the evolution of the shape and position 
of the peaks with the system's size (section l4.3p . 

4.1. The q = sector 

We have plotted in figure [10] — left the probability density at g = as a function of 
T for all our lattices. There clearly is an enveloping curve in the region T < with 
a decreasing, but positive, value of -P(O). In a mean- field setting [8] we expect this 
probability density to go to zero linearly in T. In order to check this, we have plotted 
P(0)/T against T in figure fTOl — right. As we can see, this expectation is fulfilled. For 
a similar study see [63]. We remark that the seemingly out of control value of -P(O) for 
our lowest temperature in L = 8 is an artifact of the binary nature of the couplings (a 
finite system always has a finite energy gap). Indeed, in [61], the finite size behaviour of 
P(0) for the Edwards- Anderson model with binary couplings was studied as a function 
of temperature. Finite-size effects on P(0) turned out to be stronger close to T = 
than at finite temperature. 

From a droplet model point of view, Moore et al. [21] have argued that the apparent 
lack of a vanishing limit for P(0) in numerical work in the 1990s was an artifact of critical 
fluctuations. In fact, at Tc, -P(O) diverges as L^^'^ while droplet theory predicts that, 
for very large lattices, it vanishes as L~^, with ( ~ 0.2, for all T < Tc. These authors 
rationalise the numerical findings as a crossover between these two limiting behaviours. 
However, a numerical study at very low temperatures (so the critical regime is avoided) 
found for moderate system sizes a non- vanishing P{0) [63]. Furthermore, we compute 
in section 14.31 a characteristic length for finite-size effects in the spin-glass phase, which 
turns out to be small at T = 0.703. 
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Figure 11. Infinite volume extrapolation of the Binder parameter at T = 0.703 and 
T = 0.625 and fits to the behaviour expected in the RSB, eq. (|44)) . and droplet, eq. p3|) . 
pictures. See tabled For the experimentally relevant scale of L = 110 (dotted vertical 
line, see section [6]) both fits are well above the B = 1 value of a coarsening system. 



4-2. The Binder cumulant 

We have plotted the Binder cumulant (j6]) for T = 0.625, 0.703 as a function of the system 
size in figure ITTl As discussed in section [2l2l the evolution (and thermodynamical limit) 
of this observable is different in the droplet and RSB pictures: 

Droplet : 5(T; L) = 1 + aL'^ , (43) 
RSB: B(T-L) = c + dL-^'^, (44) 

where l/z> = 0.39(5) [26]. Since it is compatible with our best estimate for the replicon 
exponent, 6'(0) = 0.38(2), we prefer to use the second, more accurate value (there is some 
analytical ground for this identification |26]). We will attempt to distinguish between 
these two behaviours by fitting our data to ( l43l) and (jH]). 

These two-parameter fits are plotted in figure [TT] and the resulting parameters are 
gathered in table |2l In the case of the RSB fit, Eq (1441) . we have included two error 
bars: the number enclosed in parentheses ( ■ ) comes from the statistical error in a fit 
fixing l/i> to ^(0) and the one inside square brackets [■] is the systematic error due to 
the uncertainty in ^(0). 

As it turns out, both fits have acceptable values of psr degree of freedom (d.o.f.). 
However, the evolution of B with L is very slow, so in order to accommodate the limit 
value of i?(L — ;> oo) = 1 consistent with the droplet picture, we have needed a very 



Table 2. Scaling of the Binder parameter and fit to the behaviour expected in the 
droplet, eq. (|43l), and RSB pictures, eq. (|44|). 



Droplet fit RSB fit 



T 


xVd.o.f. 


a 


c 


xVd.o.f. 


c 


d 


0.703 
0.625 


3.78/3 
2.00/2 


0.312(17) 
0.289(16) 


0.110(17) 
0.134(21) 


3.44/3 
2.73/2 


1.165(12)[34] 
1.128(11)[33] 


0.186(34) [03] 
0.193(28) [03] 
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small exponent {( ~ 0.12, smaller than the droplet prediction of C ~ 0.2 [H]). On the 
other hand, according to droplet theory [11], the connected spatial correlation function 
at q = qEA decays as l/r*". A direct study [26], however, indicates that these correlations 
decay as l/r°-^. 

The reader may find it disputable, from an RSB point of view, that a single power 
law should govern finite size effects. It would be rather more natural that corrections 
were of order 

' (45, 



19ML) ^ ie{q) 

It turns out, however, that 9{q) hardly depends on q (except on the neighbourhood 
of Q'ea), see [26] and Sect. El The neighbourhood of ^ea would produce a subleading 
correction of order 

In any case, see section [6], we remark that the relevant regime for comparison with 
experimental work is L ^ 110, where both the RSB and the droplet fits predict that 
B{T,L) is well above 1 (see figure [TTl) . 

4-3. The peaks of P{q), qEA, and finite size effects 

One of the features of the P{q) about which droplet and RSB agree is the fate of 
its two symmetric peaks as we approach the thermodynamical limit. These should 
grow increasingly narrow and shift their position until they eventually become two 
Dirac deltas at g = ±gEA- The actual value of gEA is notoriously difficult to 
compute [251 |65l [66] , see, however, [26] . 

Characterising the evolution of these peaks as we increase the system size is the 
goal of this section. We start by defining qEA{L) as the position of the maximum of 
P(g; L) (since the pdf is symmetric, we shall consider all overlaps to be positive in the 
remainder of this section). Thanks to the Gaussian smoothing procedure described in 
eq. dH]), this maximum is very well defined. We compute its position by fitting the peak 
to a third-order polynomial (notice that the peaks are very asymmetric). 

In order to further describe the peaks, we will also employ the half-widths cx^^-' at 
half height [P(gW) =P(gEA(^))/2]: 

aW = |gW-gEA(i^)| (46) 

where g'-"-' < q^AiL) < q^^\ 

We have plotted these parameters as a function of temperature in figure [121 On 
table [3l we can see that the width of the peaks does decrease with a power law in L, 
although very slowly. The product aP{qEA{L)) has a small dependence on L. 

We can now extrapolate q^AiL) to find the order parameter in the thermodynamical 
limit. A finite-size scaling study [26] shows that 

gEA(L,T) = g-A(T)[l + ^], A{T) = [L^iT)]'/' , (47) 

where 1/0 = 0.39(5). Yet, as discussed after eq. ( I44l) . we prefer to identify 1/0 with 
the replicon exponent, ^(0) = 0.38(2). A disagreeing reader merely needs to double 
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Figure 12. Left: q-EA{L) as a function of the temperature. We include two different 
infinite- volume extrapolations: using the replieon exponent, eq. ()47|1 and table [H and 
the one obtained from finite-size scaling arguments in the critical region, eqs. (j49p 
and ([5^ . Right: Width of the peaks of P{q), eq. (pS]). as a function of T for all our 
lattice sizes. 



Table 3. Width a = (cr(+) + cr*"' ) /2 of the peaks in P{q) and fit to a power law 



a{L) = AL^ in the range [Li„in,32]. We also include the product (jP{q^x{L)). 





T = 


0.703 


T = 0. 


805 


L 


a 


aP{qEA{L)) 


cr 


aP{qEA{L)) 


8 


0.1177(20) 


0.1784(10) 


0.1391(25) 


0.1833(10) 


12 


0.0963(21) 


0.1740(12) 


0.1165(25) 


0.1809(12) 


16 


0.0817(16) 


0.1696(11) 


0.1001(22) 


0.1756(11) 


24 


0.0735(16) 


0.1690(12) 


0.0860(19) 


0.1728(12) 


32 


0.0668(29) 


0.1631(23) 


0.0798(34) 


0.1669(22) 


-^min 


16 




16 




xVd.o.f. 


0.43/1 




1.13/1 




B 


-0.278(28) 




-0.346(30) 





the error estimate in the extrapolation of gsA- Note that one should not attempt a 
three-parameter fit to eq. (147|) . as there are too few degrees of freedom. An independent 
estimate of l/i> is required. Similar extrapolations were attempted in [17], with smaller 
system sizes (L < 16) and a lesser control over l/v. 

We present the values of q^AiL) and the result of a fit to eq. P7|) on table HI As we 
can see, the errors due to the uncertainty in the exponent, denoted by [■], are greater 
than those caused by the statistical error in the individual points, (■). In fact, our 
data admit good fits for a very wide range of values in l/z>. For instance, if we try to 
input the value of the exponent obtained in the droplet-like extrapolation of the Binder 
parameter, ( ~ 0.12 (see eq. f H3|) and table |2]), we still obtain a good fit, even though 
the extrapolated value for ^ea is almost zero at T = 0.703 and negative at T = 0.805. 
Therefore, using the droplet exponent ( the spin-glass phase would be non-existent. 

Also included in table Hlis the confidence interval for this observable computed from 
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Table 4. Extrapolation to infinite volume of q^A{L,T) using the replicon exponent, 
eq. (|47| . We also include the confidence interval previously obtained in a non- 
equilibrium study |25j . 





L 


T = 0.703 


T = 0.805 




8 


0.82461(83) 


0.7818(11) 




12 


0.79333(85) 


0.7412(11) 




16 


0.77300(75) 


0.71681(95) 




24 


0.74027(71) 


0.67905(83) 




32 


0.7174(14) 


0.6535(16) 






16 


16 


xVd.o.f. 


1.83/1 


0.98/1 




Qea 


0.538[11](6) 


0.447[12](6) 


Bounds from 


[25] 0.474 < gEA < 0.637 0.368 < qea < 0.556 






Table 5. Determination of Lc in eq. (|47|) for several temperatures below Tc. Errors 






are given as in table |4l The characteristic length Lc{T) scales as a correlation length 






when T approaches Tc {v « 2.45 from [S^). We warn the reader 


that the x^/d.o.f. for 






the fits at r = 0.85 and 0.90 are, respectively, 2.6/1 and 2.7/1. 




T 




lJ 




0.703 




1.253[10](32) 1.78[4](11) 


0.197[4](13) 


0.75 




1.448[12](34) 2.58[6](16) 


0.210[4](13) 


0.805 




1.731 [14] (44) 4.08[9](27) 


0.221[5](15) 


0.85 




2.023[16](54) 6.09[13](42) 


0.222[5](15) 


0.90 




2.514[21](66) 10.63[22](71) 


0.230[5](15) 



non-equilibrium considerations in [25] . Notice that the equihbrium values are much more 
precise, but consistent. The extrapolations included in this table (and analogous ones 
for other values of T) are plotted on figure [121 

We remark that the estimate of ^ea from eq. (H7|) is fully compatible with the results 
of a Finite-Size Scaling analysis of the conditional correlation functions |26j . 

Interestingly enough the estimate of ^ea provides a determination of the correlation- 
length in the spin glass phase. The reader might be surprised that a correlation length 
can be defined in a phase where correlations decay algebraically. Actually, finite size 
effects are ruled by a crossover length Lc{T) [67], that scales as a correlation length (i.e. 
Lc(T) oc (Tc — T)'"). In fact, one would expect qEAiT, L) /qEA{T) = 1 + h[L/Lc(T)]. 
The only thing we know about the crossover function is that it behaves for large x as 
h{x) ~ x^^/*^. Making the simplest ansatz h(x) = x~^l^ ^ the amplitude for the finite- 
size corrections in eq. (147]) can be interpreted as a power of the crossover length Lc{T), 
table [5l We note that our determination of L^iT) really scales as a bulk correlation 
length, with Tc and u taken from [32]. It turns out to be remarkably small at T = 0.703. 

The above argument tells us that good determinations of q^AiT) are possible. 
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provided that L ^ Lc{T). Yet, finite size scaling can be used as well to extrapolate 
(1ea{T,L) to the large-volume limit, even closer to where L becomes smaller 
than Lc- This somehow unconventional use of finite size scaling was started in 
Refs. [681 ESI [TOl [71], and has also been used in the spin-glass context [H [72]. Most 
of the times, these ideas are used in the paramagnetic phase, but we show below how 
to implement them in the low-temperature phase. 
Close to Tc, we know that 

q^^iT) = A(Te - Tf[l + /i(T, - rr + •••]• (48) 

We have excellent determinations of Tc and f3 from the work in [32], so we need only to 
estimate the amplitude A. In fact, Wegner's confluent corrections {T^ — T)^'^ are small 
close to Tc. To proceed, we note that finite-size scaling tells us that 

qEA{L,T) = L-^l''F{x)[l + L--G(x) + ...], x = L^'^iT, - T), (49) 

where the critical exponents are (from [32]). 

z/ = 2.45(15), /3 = 0.77(5), w = 1.0(1). (50) 

In order to connect eq. fH9l) with the infinite-volume limit in eq. pHj) the asymptotic 
behaviour of the scaling functions F{x) and G{x) must be for large x 

F{x) ~ x^, G{x) ~ x^^r (51) 

The resulting scaling plot is represented on figure [131 Varying the values of Tc and 
the critical exponents inside their error margins does not make significant changes in 
the plot. Notice how the curves collapse for small values of the scaling variable x and 
large L, but how for our lowest temperatures scaling corrections become important. In 
fact, eq. (149|) implies that when the temperature is lowered away from Tc the amplitude 
for scaling corrections grows as x'^'^ ~ x^'^^. 

In order to estimate the amplitude A we shall concentrate on the small-x region 
where finite-size scaling corrections are smallest. Disregarding scaling corrections in 
eq. (I19D, 

{qEAiL,T)L^/''f^ = F{xy/^ x. (52) 



The inset of figure [T3] shows that we reach this asymptotic behaviour for L > 24. Then, 
using the simplest parameterisation, F{x) = {X^^^x + B)^, 



qEA{L,T) = X{T,-Ty 



^ Ai//3(Tc-T)Li/- 



(53) 



We can fit our T = 32 data for x < 0.4 (where the curves for L = 24 and L = 32 
are compatible) and use the resulting value of A to extrapolate in eq. fl53p to infinite 
volume. This extrapolation is represented as a function of T on figure [121 It is clear 
that this critical extrapolation differs with the extrapolation from (1471) at most by two 
standard deviations. The difference, if any, could be explained as Wegner's confluent 
corrections. However, to make any strong claim on confluent corrections, one would 
need to estimate the error in the critical extrapolation. Unfortunately, we have found 
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Figure 13. Scaling plot of y = qEA{L, T)L^/^ in the critical region below Tc, following 
eq. (jini) and using the values given in [35] for the critical exponents and T^. Inset: 
Close-up of the region near Tc in the representation of eq. ([52]) . showing a linear 
behaviour for large L. 

that this error estimate is quite sensitive to the statistical correlation between Tc, v, 
and /9 (as far as we know, the corresponding covariance matrix has not been published). 

One could be tempted to compare eq. ( l53i) with eq. ( H7I) and conclude v = v. We 
observe that, at the numerical level, v = 2.45(15) [32] and = 2.6(3) [26]. However, we 
do not regard this as fireproof. Indeed, it is a consequence of our somewhat arbitrary 
parameterisation F{x) = {X^^^x + B)^. To investigate this issue further, the small- 
X region is not enough. One is interested in the asymptotic behaviour of F{x) for 
large x where unfortunately corrections to scaling are crucial. A careful study of the 
crossover region can be done only by considering corrections to scaling both at the 
critical temperature (at q = 0) and below the critical temperature (at q = ^ea)- 

Finally, the reader could worry about the applicability of fHSj) well below Tc. The 
issue has been considered recently within the framework of droplet theory [73]. It was 
found that (148|) is adequate for all T < Tc (actually, no Wegner's scaling corrections 
were discussed in [73]). Thus, the fact that our data are describable as scaling behavior 
with leading Wegner's correction does not imply that they are not representative of the 
low temperature phase. 

5. Conditional correlation functions 

Let us consider the conditional spatial correlation function C4(r|g), eq. ( fT4l) . A 
thorough study in the Fourier space is performed in [26] . Here, we provide some 
complementary information, concentrating on real space and considering as well the 
statistical fluctuations on the correlators. 
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Figure 14. Spatial correlation function Ca,^ {r, 0, 0)\q = O) at T = 0.703. Wc show on 
the right panel a rescaled version using the replicon exponent 6 = 0.38 and the scaling 
variable r/L. 




Figure 15. Subtracted correlation function, eq. ([54| . in units of l/L^^*^) as function 
of . We took the non-equilibrium determination of the replicon exponent, d{0) = 
0.38(2) [25]. 



We first concentrate on g = 0, the region where the droplet and RSB theory most 
differ. In figure [T3] — left we show C4^{r\q = 0) for T = 0.703, which is seen to tend 
to zero for large r. Furthermore, if we use the droplet scaling of eq. ( I20l) . we see that 
we need to rescale the correlation function by a factor L^^^\ with 9{0) = 0.38(2) the 
replicon exponent, in order to collapse the curves. 

As for other values of g, we may consider the differences 

Cir = L/A\q) - C,{r = L/2\q) ~ , (54) 

where the subtraction takes care of the large-r background in C/i{r\q). As we show in 
figure [T5| the subtracted correlation function scales in the range q^ < 0.2 as L~^^^\ This 



Nature of the spin-glass phase at experimental length scales 



32 





Figure 16. Plots of the conditional variance at fixed q of Qunk and C4(r) at T = 0.703, 
rcscalcd by appropriate powers of L (we chose exponents that provided a good scaling 
at q ~ 0). The abcissas correspond to q in units of q^j^{L,T ~ 0.703). 



implies that the connected correlation functions Ci{r\q) — q^ decay algebraically for large 
r (a similar conclusion was reached in [I9]). On the other hand, for = q^p^ 0.3, 
the exponent 9{q) is definitively larger than ^(0) (a detailed analysis indicates ^(^ea) ~ 
0.6 |26]). The crossover from the scaling ^{r = L/ A\q) - Ci{r = L/2\q) ~ to 
Ci{r = L/A\q) — Ci{r = L/2\q) ~ i/L^ii^A) Qg^^i be described by means of Finite Size 
Scaling [26] . 

Recalling that (Qiink) = C4^{r = 1), we can consider the spatial correlation as a 
sort of generalisation of the link overlap. In this sense it is worth recalling that in a 
mean-field setting fixing also fixes Qunk- In a three-dimensional RSB system one 
would, therefore, expect the conditional variance Var(Qiink|5'), eq. ( iTTl) . to tend to zero 
for large lattices [17]. The first panel of figure [T6] demonstrates that this is the case in 
our simulations, where we find that Var((5iink|Q') ~ L~^^'^. We can extend this result 
to r > 1 by considering the conditional variances of C4. Notice that, unlike Qimk, C4 
is already defined as an averaged quantity in eq. (fT3|) and not as a stochastic variable, 
so speaking of its variance is either trivial or an abuse of language. However, to avoid 
clutter, we have maintained the notation Var(C4(r)|g), as its intended meaning is clear. 
These are are plotted in figure [161 where we see that they decrease even faster than 
Var((5iink|?), with a power of L that does not seem to depend on r. 
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6. Non-equilibrium vs. equilibrium 

In reference [21], we suggested the existence of a time- length dictionary, relating results 
in the thermodynamical limit at finite time t^ with equilibrium results for finite size 
L. The matching for T = 0.7 was L ^ 3.7^(tw), where ^(tw) is the coherence length 
at time t^,. The comparison there was restricted to L < 20. The expectation value 
E{Qhnk\<l) was confronted with the correlation function C2+2{r = 1), recall the definitions 
in section 12.41 We also predicted that the equilibrium data for L = 33 would match our 
non-equilibrium results for t„ = 2^^. Using the same time- length dictionary our L = 32 
simulations would correspond to t„ ^ 2^^ and those for L = 24 would correspond to 



Now, recalling that (Qiink) is merely C4(r = 1), it is natural to extend this 
correspondence between C4(r) and C2+2{r) to r > 1. Of course, care must be exercised 
because C4 in a finite lattice cannot be computed beyond r = L/2, while C2+2 is defined 
for arbitrary r. However, the matching is very accurate, even for r dangerously close to 
L/2, see figure [T71 It is interesting to point out that the off-equilibrium results of [2l] and 
our equilibrium simulations have similar precision, even though the latter required about 
twenty times more computation time on Janus, not to mention a much more complicated 
simulation protocol. In this sense we arrive at the conclusion that simulating the 
dynamics may be the best way to obtain certain equilibrium quantities. On the other 
hand, only the equilibrium simulations give access to the crucial C(t,t^) = physics. 

We may now wonder about the experimentally relevant scale of one hour (t^ ~ 
3.6 X 10^^, taking one MC step as one picosecond [1]). Assuming a power-law behaviour, 
^{tw) = Atl!^^'^\ with z(0.64Tc) = 11.64(15) [25], we conclude that the correspondence 
is 1 hour < — y L ^ 110. Note, see for instance figure [TH that L = 110 is close enough 
to L = 32 to allow a safe extrapolation. 

Let us finally stress that the modified droplet scaling for ^(t^) [71] would predict 
that one hour of physical time would correspond to equilibrium data on L even smaller 
than 110. Indeed, according to these authors the time needed to reach some coherence 
length ^(t^) grows as 



where Tq is the microscopical time associated to the dynamics; Zc is the dynamical critical 
exponent computed at the critical point; tp is the exponent that takes the free energy 
barriers into account (from the dynamical point of view) and Y(T) = Yo{l — T/Tc)'^'^, 
with the u exponent being the static critical exponent linked to the coherence length. 
Near the critical point Y{T) — )■ and the power law critical dynamics is recovered. On 
the other hand, if we stay below T^,, Eq. fl55|) predicts an algebraic grow of t„ with C,{t„) 
only for very small coherence lengths. However, as the coherence length grows, the time 
needed to reach it diverges exponentially on ^(t„). 



t 



w 




(55) 
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Figure 17. Equilibrium C4{r\q) as a function of (lines) for L = 24 (left) and L = 32 
(right) lattices at T = 0.703. We compare with non-equilibrium data from [25] (points) 
of C2+2(f, t,tw) as a function of C'^{t,t^) for i„ ~ 2^^ (left) and t„ ~ 2^^ (right), (see 
section [2^ for definitions). The errors in both sets of data are comparable, and smaller 
than the point size. 



7. The link overlap 

We shall address here three separated problems: overlap equivalence (section I7.ip . 
replica equivalence (section W72\\ . and the scaling of the link susceptibility (section 17131) . 

7.1. Overlap equivalence 

As we have discussed previously, it has been proposed [H], [17] that attention should 
be shifted from the spin-overlap (the primary object for mean-field systems) to the 
link overlap (the would-be primary object below the upper critical dimension). Two 
requirements should be met for this change of variable to be feasible: 

(i) The conditional variance Var(Qiink|'?) must vanish in the large L limit. 

(ii) The conditional expectation ^{Qw-nklq) should be a strictly increasing function of 

The scaling with L of Var((5iink|Q'), section[5l does suggest that the first requirement 
holds. We shall investigate here the second requirement. We remark that the RSB theory 
expects it to hold, while droplet expects it not to. Furthermore, this point is actually 
the only disagreement between the RSB and the TNT picture. In fact, RSB expects the 
derivative dE((5imk|Q')/dq'^ never to vanish. On the other hand, TNT supporters expect 
this derivative to scale as L^'~^, where Ds represents the (would be) fractal dimension 
of the surface of the spin-glass domains. In D = 3, D — ^ 0.44 [16]. 

To estimate the derivative dE((5imk|9)/dg^, we observe that £'(Qiink|Q') is an 
extremely smooth function of (see the r = 1 curves in figure IT7|) . Hence we can 
attempt a polynomial fit: 

m 

E(Qiinklg) - E(gu„k|g = 0) = J2 • (56) 

k=l 
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Table 6. CoefBcients " in the fit to eq. ([56]). for various orders of the 
fitting polynomial, and T = 0.703 and 0.625. This coefficient is interpreted as 
[dE((5iink|'i')/d(j'^] ^2^o- ^® report as well the results for fits of the form C2^^ = A/L + c 
(centre) and 4"*^ = B/L^'^^ + d (bottom). For both fits, we also provide the 
extrapolation to L=: 110 which, according to the time- length dictionary, corresponds 
to the experimentally relevant length scale. 







T = 0.703 






T = 0.625 




L 


'-'2 






^2 


^2 


^2 


8 


0.403(5) 


0.405(16) 


0.43(3) 


0.414(7) 


0.423(19) 


0.45(4) 


12 


0.317(5) 


0.321(14) 


0.35(3) 


0.331(6) 


0.335(18) 


0.36(3) 


16 


0.271(4) 


0.262(11) 


0.26(2) 


0.282(6) 


0.275(16) 


0.28(3) 


24 


0.224(5) 


0.222(15) 


0.22(3) 


0.231(5) 


0.220(14) 


0.22(3) 


32 


0.199(6) 


0.201(18) 


0.20(4) 








xVd.o.f. 




0.57/3 






0.46/2 




A 




2.23(21) 






2.46(27) 




c 




0.129(16) 






0.121(21) 




L = 110 




0.149(14) 






0.143(19) 




xVd.o.f. 




2.39/3 






0.18/2 




B 




1.45(11) 






1.32(15) 




d 




-0.06(3) 






-0.11(5) 




L = 110 




0.082(20) 






0.058(28) 







Table 7. C(r = l|g) for g = and q = gEA for all our system sizes at T = 0.703. For 
each L, we include the correlation coefficient between both values of q. Specifically, for 


two quantities A and B, TZab = {{A) 


- {A))m (B))/^ 


/{{A)-{A)r m-{Bw 


L 


cm 


C(l|gEA) 


n 


8 


0.46138(82) 


0.57253(33) 


0.134 


12 


0.51649(71) 


0.60390(28) 


0.051 


16 


0.54552(60) 


0.62089(22) 


0.060 


24 


0.57573(77) 


0.63742(17) 


-0.119 


32 


0.59131(94) 


0.64579(24) 


0.063 



In particular, the coefficient provides an estimate of d-E'(Qiink|Q')/dq'^ at = 0. 

Playing with the order 2m of the polynomials, one can control systematic errors. Mind 
that it is very important to fit the difference E(Qiink|q') — E(<5iink|Q' = 0), which, due to 
statistical correlations, has much reduced statistical errors. On the other hand, data for 
different q are so strongly correlated that standard fitting techniques are inappropriate. 
We thus used the approach explained in ref. [25j. The results, see table [6l indicate that 
c^2^ offers a reasonable compromise between systematic and statistical errors. 

Once we have the derivatives in our hands, we may try to extrapolate them to 
large L by means of an RSB fit {A/ L + h, middle part of table [6]) or using a TNT fit 
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[B/L^-^^ + (i, bottom part of tabled]). The two functional forms produce a reasonable 
fit. As expected, the 1/L extrapolation to L = oo yields a non- vanishing derivative, 
while the extrapolation suggests that, for large L, E((5iink|Q') is constant as 

varies. We remark as well that the very same conclusion was reached in the analysis of 
the non-equilibrium temporal correlation functions [25] . 

However, we have far more accurate data at our disposal than the derivative 
d£'(Qimk|?)/d5^ at = 0, namely the correlation functions themselves. In table [7] 
we give our estimates for C{r = l\q = 0) and C(r = l\q = 0.523 ^ Q'ea)- According 
to a TNT picture of the SG phase, the two correlation functions should be equal. As 
the reader can check, an infinite volume extrapolation as is unbearable for both 

correlation functions (even if we discard the two smallest sizes). The same conclusions 
hold substituting L by 

£ = tt/ sin(7r/L) , (57) 

which is more natural for lattice systems. Yet, it could be argued that our data are 
preasymptotic. Hence, we may try a TNT extrapolation including scaling corrections. 

C{r = l|g) = Coo + A,L-°-^'^(l + B,L-y) . (58) 

We have performed a joint fit of the data on table [7] to eq. fl58|) . The fitting parameters 
were the four amplitudes Aq, -Bo, ^0.523 and -B0.523, the common scaling corrections 
exponent y and the common large-L extrapolation Coo- We take into account the (almost 
negligible) correlation in data for the same L by computing with the covariance 
matrix, which can be reconstructed from the data on table [71 The result is (notice the 
highly asymmetric errors) 

Coo = 0.677^°:°^^^ 2/ = 0.57lg^, xVd.o.f. = 9.1/4. (59) 

Were the functional form in eq. (158!) correct, the probability of being even larger than 
we found would be only 6%. 

On the other hand, in an RSB setting, one would expect C{r = l\q) to scale as 1/L, 
with a g-dependent infinite volume value Coo{q)- Indeed, if we fit the data on table [71 to 
C{l\q) = Coo(g) + A/i we obtain 

C^(q = 0) = 0.6349(8), xVd-o.f. = 3.63/3, (60) 
Coo(g = gEA) = 0.6711(2), xVd.o.f. = 2.86/3. (61) 

We note as well that [Coo{q = ^ea) ~ Cooil = ^)]/<1ea ~ 0.132, in fair agreement with 
the 1/L extrapolation for the derivative in table [6l 

However, more important than the extrapolation to L = 00 is the extrapolation to 
L = 110, the length scale that, for T = 0.7, matches the experimental time scales. For 
T = 0.625, L = 110 is surely larger than the relevant length scale but, unfortunately, 
the time-length dictionary at such a low temperature still needs to be tuned. As it 
can be seen in the middle and bottom parts of table [6l the two extrapolations yield a 
non- vanishing derivative. 
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Thus, whichever the standpoint adopted, the conclusion is identical for RSB and 
TNT theories: at the experimentally relevant length scales, overlap equivalence can be 
assumed. 

7.2. Replica equivalence 
We consider now the ratio 



-Klink — ^ 2 ' K^'^J 

(Qvmk) ~ (Qlink) 

defined in section 12.51 As was explained there, the RSB theory expects it to reach 
a constant value 2/3 below Tc, whereas the droplet and TNT theories lack a definite 
prediction. Our numerical data fit very well the RSB expectation (see figure [T8]-left). 

Besides, we can also study a similar ratio, in which the mean-field substitution 
Qiink is performed: 



Rq^ = = =2 ■ (63) 

Overlap equivalence suggests that Rg2 approaches 2/3 in the large L limit (again 
neither the droplet nor the TNT theories have a definite prediction). Our data at 
low temperatures seem compatible with the 2/3 expectation, see figure [T8]-right. On 
the other hand, the convergence to the thermodynamic limit seems fairly slower close to 
Tc. We recall that a previous computation concluded as well that violations of Rq2 =2/3 
are due to critical fiuctuations [9]. 

7.3. Link susceptibility 

We show in figure [TOl-left the link susceptibility Xunk, eq. (j27l) . as a function of 
temperature for different lattice sizes. It is clear enough that this susceptibility is 
divergent in the spin-glass phase and that the lower the temperature, the more violent 




0.4 0.8 1.2 1.6 ■ 0.4 0.8 1.2 1.6 

T T 



Figure 18. The ratios i?iink, eq. (PT|) . (left panel) and i?q2, eq. (right panel) 

versus T for the different system sizes. The replica equivalence property implies that, 
in an RSB system below Tc, i?"""^ = 2/3 in the large-L limit. Recall that Tc 1.1. 
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Figure 19. (Left) Susceptibility xunk vs. temperature for the different system sizes. 
(Right) Behaviour of xunk with L for different temperatures. Lines are power-law fits. 
The effective exponents found in each of these fits is reported in the legends. 





Table 8. 


Ratio S'ji^™^ , eq. (p^ , for all our 


lattice sizes at T = 0.625, 0.703, using the 




coefficients from table |6l 








T 


= 0.703 


T = 


0.625 


L 


c(2) 
^link 


^link 


o(2) 
^link 


0(4) 
'-'link 


8 


0.838(21) 


0.846(67) 


0.859(29) 


0.897(81) 


12 


0.777(25) 


0.797(70) 


0.801(29) 


0.821(88) 


16 


0.755(22) 


0.706(59) 


0.766(33) 


0.729(85) 


24 


0.776(35) 


0.76(10) 


0.745(32) 


0.675(86) 


32 


0.816(49) 


0.83(15) 







the divergence. Hence, it is clear that this particular effect is not due to critical 
ffuct nations. 

We perform a more quantitative study in figure [T9]-right. As discussed in section |23| 
according to RSB theory, one would expect xunk ~ in the SG phase. 

We find evidence of a critical divergence. At and above T^, our data grow very softly 
with L (at T = 1.3 ~ 1.17Tc, data seem to reach a limiting value). However, below Tc, 
we observe an effective exponent that grows when we lower the temperature We observe 
that the effective exponent, for our lattice sizes and temperatures, has already grown 
beyond the Chayes bound of D/2 but still has not reached the RSB expectation of D. 
Note that no existing theory of the spin-glass phase can accommodate a temperature- 
dependent exponent. Therefore, the most economic scenario is that our lattice sizes are 
not large enough, so we are still in a preasymptotic regime for this quantity. 

Let us take a slightly different point of view. Rigorous theorems discussed in 
section [2751 tell us that, eq. (l30l) . if 

lim xiin^/L'' > , (64) 
also the width <Jq^.^^. of the probability density function for Qunk, 

^Q,,, = (QLJ - W^' , (65) 
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will be non- vanishing in the thermodynamic limit (see also section I7.2p . It is very 
important that the converse statement also holds. 

Now, using the identity ( lT2l) . we can split up this variance in two different 
contributions: 



4un.= J dg P(g) (^Var(Qunkk) + [E (gunkk) - (giink)J j ■ (66) 

Since Var((5imk|?) scales as L"-^/^ (see [17] and figure [T6|) . only the second term may 
survive the large L limit. 

This suggests the definition of a modified link susceptibility: 

2 



_ JZ, dq P{q) [E (Qnnklg) - (giink) 

Xlink „ 



_ (67) 
According to RSB theory, ximk should scale as whereas it would not diverge as 



violently in a droplet or TNT scenario. The rationale for dividing out the (g^) — (g^) 
can be found in eq. flSB|) . Assuming that the lowest order polynomial is adequate, one 
finds that (of course, the particular value of the index m should be immaterial) 

Xnnk^L^ict'^f. (68) 

Hence, the TNT theory would expect Xnnk/L^ to tend to zero, just because it predicts 
that in the large-L limit c^"^'' = 0. Note that the droplet theory would predict a 



vanishing xvmk/ for a different reason, namely because they expect that {q^) — (g^) 
should vanish. 

Let us check to what extent the estimate (l68l) is accurate. We show on table [8] the 
ratios 

(2) _ L [C2 J (4) _ L [C2 J 

•^link - ' '^link - • l^JJ 

Alink Alink 



Referring again to eq. (l56l) . it is clear that the contribution linear in g^ explains a large 
fraction of Ximk, and that this fraction is not likely to vanish in the large-L limit. 

Hence the question of whether ximk diverges as or not, turns out to be strictly 
equivalent to that of overlap equivalence that we discussed at length in Sect. 17.11 Our 
interpretation is that the effective scaling in figure ITOl-right is mostly due to strong finite 
size effects in c^"^\ Under this light, the effective exponents reported in figure flQl-right 
are preasymptotic. In fact, the ratio A/cis large {c^"^\l) = c+A/L, see table[6]), which 
tells us that for xunk and related quantities finite volume corrections are particularly 
large and naive power law fits may give wrong results. 

Let us conclude this section by checking how these quantities behave in a 2D Ising 
ferromagnet (i.e. with no disorder built in). Although this model is clearly too simple, 
it is also true that, up to our knowledge, the quantities investigated here have not 
been looked at before. Hence, it is interesting to see what happens even in this simple 
case. We use two replicas to compute Xunk- Results for ximk are presented on table [9] 
for two different temperatures below the critical temperature Tc. There we can see 
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that xiink approaches a hmiting value when L grows. Furthermore, the Umiting 

value decreases when lowering the temperature away from T^. Hence, a divergent link 
susceptibility below Tc is something that should not be taken for granted. 





Table 9. xiink in the 2D-Ising model (Tc = 2/log(l + ^2) s 


i 2.26918531...). 


L 


T = 0.992rc 


T = 0.986Tc 


8 


7.16(1) 


7.091(8) 


12 


8.776(6) 


8.614(15) 


16 


10.247(12) 


9.868(9) 


24 


11.47(2) 


10.51(2) 


32 


12.058(15) 


10.639(12) 



8. Conclusions 

We have obtained equilibrium configurations of the Ising spin glass (-0 = 3, ±1 Edwards- 
Anderson model) on large lattices at low temperatures (T = 0.64Tc for L = 32, 
T = 0.56Tc for L = 24, and even lower temperatures for smaller systems, see tabled]). 
This unprecedented computation has been made possible by the Janus computer. 
However, the parallel tempering had never before been put to such stress, and we have 
devoted a large effort to convince ourselves that thermalisation was achieved. New 
thermalisation tests were devised. Furthermore, a new simulation strategy had to be 
employed: the simulation time needs to be tailored sample by sample (for one cannot 
afford adopting worst-case parameters). 

The main conclusion we draw is that the correspondence between equilibrium 
results and non-equilibrium dynamics (much easier to compare with experimental work), 
is deeper than anticipated. In fact, one can construct a time- length dictionary, such that 
equilibrium correlation functions on finite systems match non-equilibrium correlators at 
finite time (but infinite system size). The evidence for this correspondence consists 
of: (i) quantitative comparison of the spatial correlation functions and (ii) the analysis 
of overlap equivalence on equilibrium (this work) and non-equilibrium settings |25]. In 
addition, there is a remarkable coincidence between the replicon exponent obtained from 
equilibrium methods [26], and from non-equilibrium dynamics [2^ [25]. 

The unavoidable consequence of this time-length correspondence is that the system 
size that is relevant for the experimental work (time scales of one hour, say) at T = 0.64Tc 
is not infinite, but L = 110. Note that this correspondence was obtained assuming a 
power-law growth with time of the spin-glass coherence length in experimental samples. 
Should the modified droplet scaling for ^(t^) hold [71], the relevant equilibrium system 
size would be even smaller. It is obvious that extrapolating numerical data from L = 32 
to L = 110 is far less demanding than extrapolating them to infinite size. All such 
extrapolations in this work (even those assuming droplet scaling) were conclusive. The 
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only effective theory that is relevant at experimental time scales is Replica Symmetry 
Breaking. 

However, the question of whether RSB is only an effective theory in = 3 or 
a fundamental one does not lack theoretical interest. We have attempted several 
extrapolations to infinite system size in this work, finding that droplet theory is ruled 
out, unless a change of regime arises for system sizes much larger than our reached 
L = 32. We remark that in Sect. 14.31 we have numerically determined a crossover 
length that rules finite size effects. As expected for a large enough system, it scales 
with temperature as a bulk correlation length. However, on the basis of numerical data 
alone, one can never discard that new behaviour might appear for much larger system 
sizes, irrelevant for current experimental work. 

We found three contradictions with droplet theory. First, in order to have a trivial 
Binder cumulant, finite-size corrections had to be of order ~ L~°'^^. Such finite size 
corrections would imply a vanishing, or even negative, spin- glass order parameter ^ea- 
Second, according to droplet theory (see [H], page 139) finite size corrections ~ L~°-^^ 
imply that the connected spatial correlation function at q = Qea decays as 1/r^'^^. A 
direct estimate indicates that, at g = gEA, correlations decay as l/r^-^ [26]. Third, the 
probability density function P{q = 0) does not decrease with increasing system size (a 
similar conclusion was reached in [631 ES] ) • 

Our analysis of overlap equivalence is compatible with the RSB picture, without 
invoking sophisticated finite-size effects. On the other hand, the statistical likelihood 
for TNT theory, as formulated in [16], has been quantified to be 6%. In any case, 
TNT scaling predicts that for L = 110 the surface-to- volume ratio of the magnetic 
domains is still of order one (in agreement with RSB). In addition, we find that replica 
equivalence is consistent with the RSB picture (while TNT lacks a definite prediction). 
Furthermore, the link susceptibility, Xunk? is definitively divergent in the spin-glass phase 
(since the divergence is stronger the lower the temperature, its origin is obviously non- 
critical). We are aware of no argument in TNT theory implying the divergence of the link 
susceptibility. On the other hand, RSB theory does require a divergent Xunk- However, 
RSB demands a scaling Ximk ~ L^- Such growth regime has still not been reached for 
our system sizes, although we have identified the origin of this preasymptotic behaviour. 

A final lesson from the present numerical study is that careful non-equilibrium 
simulations [25] are almost as rewarding as the equilibrium work. Indeed, our previous 
non-equilibrium study [2^ [25] reached a time scale that corresponds to the present 
equilibrium L = 32 simulation. Yet, the numerical effort to obtain the data in Fig. [17] 
has been larger by, roughly, a factor of 20 in the case of the equilibrium work. It is true 
that the equilibrium approach allows to investigate directly the crucial q = region, 
where in the nonequilibrium case one would need to rely on difficult extrapolations to 
infinite time. However, we do not think that there is much road ahead for equilibrium 
studies, due to the failure of the parallel tempering algorithm. Indeed, see table [H it 
takes about 3.5 times more numerical work to equilibrate 1000 samples of L = 32 at 
T = 0.64Tc than 4000 samples of L = 24 down to T = 0.56Tc. Clearly enough, the 
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temperature window accesible with the parallel tempering algorithm decreases very fast 
as the system size grows. We believe this failure to be due to a genuine temperature- 
chaos effect. However, in order to analize quantitatively the effect one needs to correlate 
the (sample dependent) temperature bottlenecks, see Fig. [l]-left, with the spin overlap 
at different temperatures. This analysis is left for future work [62] . 
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Appendix A. Our thermalisation protocol 

We have followed a three-step procedure to thermalise each sample: 

(i) We simulate for a fixed minimum length of MCS, chosen to be enough to 
thermalise most of the samples. Notice that most published parallel-tempering 
simulations stop here, assessing the thermalisation only through the time evolution 
of disorder-averaged observables. 

(ii) We discard the first sixth of the measurements and compute the integrated 
autocorrelation time, choosing the self-consistent window W of (1371) so that W > 
6rint. Using this first estimate of the integrated time, we enlarge the simulation 
until A^'hb > 22rint, always discarding its first sixth. A criterion based on Tint was 
first used in [27] . 

(iii) Now that we have a reasonably dimensioned simulation, we can compute the 
exponential autocorrelation time (which is typically bigger than, but of the same 
order of magnitude of, the integrated time). We demand that Ahb be larger than 

-■- ^ ' exp • 

This last step is the main innovation of these simulations. Notice that we have to 
perform non-linear fits in some 10^ autocorrelation functions, with a sample-dependent 
fitting range. This is a somewhat delicate procedure, so we have taken great care to 
ensure it is failsafe. 

We start by assuming that the correlation function (l38l) can be approximated by 
the sum of two exponentials: 

C{t) ~ A^e-'/^^ + A2e-*/^^ n = re.p > r^. (A.l) 

Usually, one chooses the range for such a fit manually but this is not practical here, due 
to the sheer number of correlation functions we have to study. Hence, assuming that 
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the exponential time is not much larger than the integrated one, we have used the latter 
in order to define our fitting range (notice that if A2 = 0, Ai = 1 and ti = Tint)- Our 
fitting procedure has three steps 

(a) We perform a first fit to a single exponential in the range [2rint, 3rint], from which 
we obtain an amplitude A and a time r. 

(b) Using Ti = T, T2 = t/10, Ai = A and A2 = 1 — A as a starting point we perform 
the non- linear fit to ( lA.ll) with a Levenberg-Marquardt scheme [77]. The fitting 
range is chosen as [rjnt/lO, lOrjnt]. 

(c) Sometimes T2 is very small and C(t) is indistinguishable from a single exponential 
in [rint/10, lOrint]- In these occasions the fit in step 2 fails, which can be detected 
in a number of ways (very large or even negative values for one of the Ai, absurdly 
large values for ti or even a complete breakdown of the iterative method). For these 
samples, a third fit to a single exponential is performed in the range [5rint, lOrjnt]. 
There is one exception: when one of the A^ is negative, there appears a very 
pronounced downwards fiuctuation in C{t) for large times, which can lead to 
an underestimation of r^xp- In these occasions, the third fit is performed in 

[2.5rint, STint]- 

This automatic and fully quantitative procedure works for most samples, but there are 
some potential pitfalls which may lead to our underestimating r^xp- Sometimes, the 
exponential time is much larger than Tint- This can result in a failure of the automatic 
method for two reasons: (1) as Toxp ^ Tint, the fitting ranges are no longer well adjusted 
(2) a very large Texp/rint implies a very low value for Ai. We address this problem by 
enlarging the measurement bins by a factor of 10 in case Texp > lOrint. This way, both 
Tint and Ai grow, and the fit works much better. 

The possibility also exists that C{t) may be misleading, because the simulation is 
so much shorter than the exponential time that some of the configurations have not yet 
explored the relevant minima of the free energy (i.e., the C{t) we are measuring is not 
yet the equilibrium one). This happens when some of the 4Nt configurations have not 
crossed the critical temperature in the parallel-tempering dynamics. The assumption 
here, key to the parallel tempering method, is that once a configuration spends a few 
MCS at high temperatures it becomes completely decorrelated (remember that, due 
to Janus' special characteristics, the interval between measurements is very large). To 
prevent this from happening, we measure the time thot that each configuration spends 
at temperatures greater than T^. In case any of the AN^ configurations has a value of 
thot smaller than one third of the median, we consider that the simulation is far too 
short for us to measure Texp and we simply double A^hb- Notice that this last criterion 
is unlike the others in that it is not completely quantitative. It simply detects that our 
starting point is very badly dimensioned. 

As a final test, we have increased Amin by a factor of 10 for the first 1% of the 
samples in all lattices. None of the Texp estimates changed within errors. 
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Appendix B. Unbiased estimators of non-linear functions 

Non linear functions of thermal mean values, that are computed sample by sample and 
afterwards averaged over disorder, are prone to suffer systematic errors larger than the 
statistical ones. General cures for this problem are known [611 132]. our case, the only 
such quantity is Ximk, defined in eq. ( 1271) . Since we have 4 replicas, the bias problem 
could be avoided for Ximk, see eq. (1B.6P below. However, if one decides instead to face 
it, a nice test for the statistical quality of the data is obtained. 

Indeed, consider eq. ([27]), and let [Qunk] be our Monte Carlo estimate of (Qunk) as 
computed from measurements for a given sample. The expectation value ([Qiink]^) 
is not (Qiink)^- To quantify the effect, we need some notation [511 [59]. The normalised 
equilibrium autocorrelation function for Qunk? at a given temperature and for a given 
sample, is 

_ ((Qltt*^ - (Qlink))(Q;i - (gunk))) 
Wlink/ ~ Wlink/ 

where Q[iI]^ stands for the value taken by Qimk at time s. Two characteristic time scales 
are relevant to us: 

^ t=+oo s-^t=+oo |j| A /.\ 

— - r (f^ n- — ^^=-00 I'^l'-^Qlinkl'^J ,-r. r,\ 

'int.Qlink 9 / ^ '-'Qiinkv''/' 'avg,Qiink Y^t = + oo A /,\ ■ \^-'^) 

Then a straightforward computation shows that 

l2\ _ /A^ \2 , ^Tint.QiinkKQunk) (Qlink)^] ''"avg.Qunk 



12\ //n \2 I "'iut,^!iinkL\^link/ \-«miK/ j ' avg.v^iink \ m q\ 
([C^linkJ ) = (<^link) H [I J , {0.6) 

up to corrections of order 0(e~^/^'=''p) (the exponential autocorrelation time Tcxp was 
discussed in section 13731) . This computation is performed in textbooks [HH |59] only to 
order and with a rather different aim: it provides an estimate of the (squared) 

statistical error in the Monte Carlo estimation of (Qiink)- Out interest in eq. flB.3P is 
different. It tells us that, when taking [Qiink]^ as (Qunk)^ we are incurring in bias not 
only of order but also of order l/N"^. 

It is easy to obtain bias-corrected estimators [61]: one divides the Monte Carlo 
history in two halves, four quarters, and eight eighths. Recall that we will be dropping 
in the analysis the full first half of the Monte Carlo history. Then, one computes [Qiink]2/2 
from the last half of the data as well as [Qiink]3/4 and [Qiink]4/4 from the third and fourth 
quarters respectively. Similarly, we compute [Qunkls/g^ [Qimk]l/s, [Q\mk]ys [Qiinklg/g- 
Then, eq. (lB.3p . the thermal expectation value of 

„(2) i2 [Qlink]3/4 + [Qlink]4/4 

Vlink.lincar = 2[Qlinkj2/2 ^ > (B.4) 

is (Qiink)^, up to a bias of order Tint,Qiink''~avg,Q,i„k/^^- We can do it even better: 
n(2) =-\0 

^link. quadratic o [ t^I 



[Qlink]3/4 +[Q-- ■ 

2 

1 [QlinkJs/g + [Qlink]g/8 + [Qlink]7/8 + [Qlinkjg/g 



(2) 12 [VlmkJ3/4 -r [Vlmkj4/4 

'link.quadratic ~ 2[Vlmkj2/2 ~ ^ - 



2 I \n 12 I \n 12 , \n 12 

(B.5) 
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has thermal expectation value (Qunk)^; up to corrections of order 0(e~^/'^™p). 

The fact that we have 4 real replicas offers us an alternative way of overcoming this 
problem. In fact, denoting by Q^jn^ the link overlap computed from replicas i and j, we 
have 



r^(12)^(34) ^(13)^(24) ^(14)^(23)n 
)(2) _ [VlinkVlink ~r VlinkVlink ' VlinkVli nkJ 
'link,4R ~ g 



/OV^^ link ^ link ' ^ link ^ link ' ^ link ^ link J /id c\ 

Vlink.4R ~ o ' \^-^) 



with (Qiink4R) = (Qiink)^ (wc average over the three equivalent replica pairings to reduce 
statistical errors). The comparison of the two procedures offers an interesting test on 
the statistical quality of our data, because eq.( ]B.3l) holds only for ^ '''int.Qiink' '^cxp • 
From the different estimators for the (Qimk)^ we finally obtain four estimators of 

Xiink: 



Xlink,2R =^([^?Lk]-[Qlink]2), (B.7) 



XhSr — ^([Qunk] Qlinkjinear) ; (B-8) 



Xlink,2R ,quadratic ) , (B.9) 



Xlink,4R =^([QLk]-Qi?nU)- (^-10) 
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